﻿#region Copyright

//  The author of this software is M.G. Johnson.
//  Permission to use, copy, modify, and distribute this software
//  for any purpose without fee is hereby granted, provided that
//  this entire notice is included in all copies of any software
//  which is or includes a copy or modification of this software
//  and in all copies of the supporting documentation for such
//  software.  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT
//  ANY EXPRESS OR IMPLIED WARRANTY.  IN PARTICULAR, NEITHER THE
//  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY
//  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS
//  FITNESS FOR ANY PARTICULAR PURPOSE.

#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.Optimization
{
    /// <summary>
    /// This class represents an multidimensional optimizer, which use the Hooke and Jeeves algorithm. The algorithm based 
    /// one an implementation of M.G. Johnson. For further information see http://netlib.org/opt/hooke.c.
    /// </summary>
    [Serializable]
    public class HookeJeevesOptimizer : AbstractMultidimensionalOptimizer
    {
        /// <summary>
        /// Describes if the calculation has a precision error.
        /// </summary>
        private bool precisionError;

        /// <summary>
        /// This field holds the number of needed iterations.
        /// </summary>
        private int neededIterations;

        /// <summary>
        /// Initializes a new instance of the <see cref="HookeJeevesOptimizer"/> class.
        /// </summary>
        /// <param name="function">The base function to optimize.</param>
        public HookeJeevesOptimizer(IHardMultivariateRealFunction function)
            : base(function)
        {
            this.precisionError = false;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="HookeJeevesOptimizer"/> class.
        /// </summary>
        /// <param name="function">The base function to optimize.</param>
        public HookeJeevesOptimizer(HardMultivariateRealFunction function)
            : base(function)
        {
            this.precisionError = false;
        }

        /// <summary>
        /// Gets a value indicating whether a precision error during the calculation.
        /// </summary>
        /// <value>True if the calculation has a precision error otherwise, false.</value>
        public bool PrecisionError
        {
            get { return this.precisionError; }
        }

        /// <summary>
        /// Gets the number of needed iterations to find the minimum of the specified function.
        /// </summary>
        /// <value>The number of needed iterations to find the minimum of the specified function.</value>
        public int NeededIterations
        {
            get { return this.neededIterations; }
        }

        /// <summary>
        /// Finds the minimum of the specified function by using the Hooke and Jeeves algorithm.
        /// </summary>
        /// <param name="startValues">The start vector of the finding process.</param>
        /// <returns>The minimum of the function.</returns>
        public GeneralVector FindMinimum(GeneralVector startValues)
        {
            return this.FindMinimum(startValues, 1e-5, 0.5, 1000);
        }

        /// <summary>
        /// Finds the minimum of the specified function by using the Hooke and Jeeves algorithm.
        /// </summary>
        /// <param name="startValues">The start vector of the finding process.</param>
        /// <param name="iterations">The maximum number of interations.</param>
        /// <returns>The minimum of the function.</returns>
        public GeneralVector FindMinimum(GeneralVector startValues, int iterations)
        {
            return this.FindMinimum(startValues, 1e-5, 0.5, iterations);
        }

        /// <summary>
        /// Finds the minimum of the specified function by using the Hooke and Jeeves algorithm.
        /// </summary>
        /// <param name="startValues">The start vector of the finding process.</param>
        /// <param name="precision">The precision of the result.</param>
        /// <param name="iterations">The maximum number of interations.</param>
        /// <returns>The minimum of the function.</returns>
        public GeneralVector FindMinimum(GeneralVector startValues, double precision, int iterations)
        {
            return this.FindMinimum(startValues, precision, 0.5, iterations);
        }

        /// <summary>
        /// Finds the minimum of the specified function by using the Hooke and Jeeves algorithm.
        /// </summary>
        /// <param name="startValues">The start vector of the finding process.</param>
        /// <param name="precision">The precision of the result.</param>
        /// <param name="stepSize">The start step size of the algorithm.</param>
        /// <param name="iterations">The maximum number of interations.</param>
        /// <returns>The minimum of the function.</returns>
        public GeneralVector FindMinimum(GeneralVector startValues, double precision, double stepSize, int iterations)
        {
            GeneralVector entpt = new GeneralVector(startValues.Count);
           
            this.neededIterations= this.InternalFindMinimum(startValues.Count, startValues, entpt, stepSize, precision, iterations);
            this.precisionError = (neededIterations == iterations);

            return entpt;
        }

        private int InternalFindMinimum(int nvars, GeneralVector startpt, GeneralVector endpt, double rho,
                                        double epsilon, int itermax)
        {
            int i;
            double[] delta = new double[nvars];
            GeneralVector newx = new GeneralVector(nvars);
            GeneralVector xbefore = new GeneralVector(nvars);

            for (i = 0; i < nvars; i++)
            {
                newx[i] = xbefore[i] = startpt[i];
                delta[i] = Math.Abs(startpt[i]*rho);

                if (delta[i] == 0.0)
                {
                    delta[i] = rho;
                }
            }

            int iadj = 0;
            double steplength = rho;
            int iters = 0;
            double fbefore = this.Function.SolveAt(newx);

            while ((iters < itermax) && (steplength > epsilon))
            {
                iters++;
                iadj++;

                for (i = 0; i < nvars; i++)
                {
                    newx[i] = xbefore[i];
                }

                double newf = this.BestNearby(delta, newx, fbefore, nvars);

                int keep = 1;

                while ((newf < fbefore) && (keep == 1))
                {
                    iadj = 0;

                    for (i = 0; i < nvars; i++)
                    {
                        if (newx[i] <= xbefore[i])
                        {
                            delta[i] = 0.0 - Math.Abs(delta[i]);
                        }
                        else
                        {
                            delta[i] = Math.Abs(delta[i]);
                        }

                        double tmp = xbefore[i];

                        xbefore[i] = newx[i];
                        newx[i] = newx[i] + newx[i] - tmp;
                    }

                    fbefore = newf;
                    newf = this.BestNearby(delta, newx, fbefore, nvars);

                    if (newf >= fbefore)
                    {
                        break;
                    }

                    keep = 0;

                    for (i = 0; i < nvars; i++)
                    {
                        keep = 1;

                        if (Math.Abs(newx[i] - xbefore[i]) > (0.5*Math.Abs(delta[i])))
                        {
                            break;
                        }
                        else
                        {
                            keep = 0;
                        }
                    }
                }

                if ((steplength >= epsilon) && (newf >= fbefore))
                {
                    steplength = steplength*rho;

                    for (i = 0; i < nvars; i++)
                    {
                        delta[i] *= rho;
                    }
                }
            }

            for (i = 0; i < nvars; i++)
            {
                endpt[i] = xbefore[i];
            }

            return (iters);
        }

        private double BestNearby(double[] delta, GeneralVector point, double prevbest, int nvars)
        {
            int i;
            double minf = prevbest;
            GeneralVector z = new GeneralVector(nvars);

            for (i = 0; i < nvars; i++)
            {
                z[i] = point[i];
            }

            for (i = 0; i < nvars; i++)
            {
                z[i] = point[i] + delta[i];

                double ftmp = this.Function.SolveAt(z);

                if (ftmp < minf)
                {
                    minf = ftmp;
                }
                else
                {
                    delta[i] = 0.0 - delta[i];
                    z[i] = point[i] + delta[i];
                    ftmp = this.Function.SolveAt(z);

                    if (ftmp < minf)
                    {
                        minf = ftmp;
                    }
                    else
                    {
                        z[i] = point[i];
                    }
                }
            }

            for (i = 0; i < nvars; i++)
            {
                point[i] = z[i];
            }

            return (minf);
        }
    }
}